Phase diagram of nuclear "pasta" and its uncertainties in supernova cores 
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We examine the model dependence of the phase diagram of inhomogeneous nulcear matter in 
supernova cores using the quantum molecular dynamics (QMD). Inhomogeneous matter includes 
crystallized matter with nonspherical nuclei - "pasta" phases - and the liquid-gas phase separating 
nuclear matter. Major differences between the phase diagrams of the QMD models can be explained 
by the energy of pure neutron matter at low densities and the saturation density of asymmetric 
nuclear matter. We show the density dependence of the symmetry energy is also useful to understand 
uncertainties of the phase diagram. We point out that, for typical nuclear models, the mass fraction 
of the pasta phases in the later stage of the collapsing cores is higher than 10-20 %. 

PACS numbers: 26.50.+x,21.65+f,97.60.Bw 
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I. INTRODUCTION 

Terrestrial matter basically consists of spherical nu- 
clei. However, such an ordinary picture, that is, "nuclei 
are spherical," might not be true in collapsing supernova 
cores just before bounce and in the deepest region of inner 
crusts of neutron stars. In dense matter close to the nor- 
mal nuclear density po ~ 0.16fm~ 3 , nuclei would adopt 
nonspherical shapes including rod and slab. Phases with 
these exotic nuclei are called "pasta" phases [l|, S| • 

The pasta phases attract the attention of many re- 
searchers in the fields of nuclear physics [3| and astro- 
physics Pasta phases have important astrophysi- 
cal effects on, e.g., neutrino opacity in supernova cores 
0; H| j neutrino emissivity in neutron star cooling 0-[I3 , 
etc. Moreover, it has been shown in our previous work [8[ 
that the pasta phases would occupy 10-20 % of the mass 
of collapsing stellar core. In such a case, the pasta phases 
could have a remarkable impact on neutrino transport in 
the core and hence success of supernova explosion. 

Equilibrium states of the pasta phases have been in- 
vestigated in many earlier works (e.g., Refs. $L Il3l - fl9l ]). 
These works have confirmed that, with increasing den- 
sity, nuclear shape basically changes in the sequence 
sphere, rod, slab, rod-like bubbles, spherical bubbles, 
and, finally, uniform nuclear matter (in some nuclear 
models, however, all of the above pasta phases do not 
appear @, SS S3)- Although these earlier works 

have studied the phase diagram of the pasta phases us- 
ing various nuclear models, the following two points are 
worth consideration. First, in these works (except for 
Ref. [13[) authors consider the above-mentioned specific 
nuclear structures and determine the equilibrium state by 
comparing the free energy among them. It is hardly pos- 



sible to know in advance whether these assumed phases 
include the true equilibrium state or there are other more 
stable states. Thus, we have to examine how the phase 
diagram is changed by relaxing this assumption. Second, 
in collapsing supernova cores, where the pasta phases 
would appear, temperature reaches typically a few MeV. 
Thermal fluctuations on the nucleon distribution are not 
completely negligible considering that the nucleon Fermi 
energy and the nuclear binding energy are from several to 
tens MeV. However, thermal fluctuations cannot be prop- 
erly incorporated by the framework employed in the ear- 
lier works such as a liquid-drop model and the Thomas- 
Fermi approximation. 



To overcome the above problems, we use quantum 
molecular dynamics (QMD) 22-25]. In these works, 
we have confirmed the pasta phases appear at zero and 
nonzero temperatures and the sequence of nuclear struc- 
tures with increasing density is the same as that in the 
earlier works. There we have also obtained spongelike 
"intermediate" phases. However, we have studied the 
pasta phases using only one specific nuclear force. We 
note that uncertainties of nuclear force, especially those 
of the surface energy and the symmetry energy, have 
larg e effects on the phase diagram at subnuclear densities 
[I6HI8I . |26|. Thus, in the present work, we shall reveal 
the influence on the phase diagram by uncertainties of 
nuclear force in the framework of QMD. 



In the followings, we set the Boltzmann constant kg 
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II. FRAMEWORK OF QUANTUM 
MOLECULAR DYNAMICS 



Models 



In our previous studies |22h24 l27j . we used nuclear 
force developed by Maruyama et at. (Model 1) [28| with 
medium equation-of-state (EOS) parameter set. In the 
present work, we also use another model by Chikazumi 
et al. (Model 2) [29] to investigate the model dependence 
of phase diagram. The Hamiltonian of both the models 
is written as, 



where K is the kinetic energy; Vp a uii is the Pauli po- 
tential, which is introduced to reproduce effects of the 
Pauli exclusion principle; Vskyrme is the Skyrme-type in- 
teractions; V sym is the symmetry energy; V^urfaco is the 
potential dependent on the density gradient; Vmd is the 
momentum-dependent potential in the form of the ex- 
change term of the Yukawa interaction; and Vcouiomb is 
the Coulomb potential. Each term is expressed as follows 



H=K + Vp au li 

+ Vskyrme + V sym + Surface + VlvlD + Vcouiomb, 
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where R^ and Pi are the centers of position and momen- and the single- nucleon densities pi(r) and pi(r) are given 
turn of the wave packet of zth nucleon and rtii, Ui, and by 
Cj (cj=l for protons and ^=0 for neutrons) denote the 
mass, the spin, and the electric charge (in units of e) 
of ith nucleon. Here p^j means the overlap between the 
densities of ith and jth nucleons, 
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with the normal width L w and the modified width L w of 
the wave packet, 



Pij = / drpi(r)pj(r), 



(9) 
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TABLE I: Parameter sets for model 1 28] and model 2 Y2& 





Model 1 


Model Z 


Cp (MeV) 


207 


115 


po (MeV/c) 


120 


120 


go (fm) 


1.644 


2.5 


a (MeV) 


-92.86 


— 121.9 


ft (MeV) 


169.28 


197.3 


T 


1.33333 


1.33333 


C s (MeV) 


25.0 


25.0 


Vsf (MeV) 





20.68 


C y J (MeV) 


-258.54 


—258.54 


Cix' (MeV) 


375.6 


375.6 


in (fm" 1 ) 


2.35 


2.35 


M2 (fm" 1 ) 


0.4 


0.4 


Ll (fm 2 ) 


2.1 


1.95 


po (fm~ 3 ) 


0.165 


0.168 



(The squared widths L\ and L\ correspond to L and L, 
respectively, in the notation of Refs. p8l[29|.) 

Parameters for the models are shown in Table HI Note 
that Vsf = for model 1, i.e., model 1 does not include 
Kiurfacc- Model parameters qo, po, and Cp in the Pauli 
potential are determined by fitting the kinetic energy of 
free Fermi gas at zero temperature. The other model 
parameters are determined to reproduce the saturation 
properties of symmetric nuclear matter [i.e., the satura- 
tion density (~ 0.16 fm -3 ), saturation energy (—16 MeV 
per baryon), and incompressibility (280 MeV)], and the 
binding energy and rms radius of the ground state of sta- 
ble nuclei. Especially these properties of heavy nuclei are 
better reproduced by model 2 than by model 1 due to the 
term V^ ur f ace [3l|. Note that V^ ur f ace is just a potential 
depending on the density gradient and is different from 
the surface energy. The surface energy comes from an 
energy loss due to the deficiency of nucleons interacting 
with each other in the region of the nuclear surface. 



B. Equations of motion 

We show equations of motion of QMD, which we 
employ to simulate the equilibrium states at zero and 
nonzero temperatures. The Hamiltonian form of the 
QMD equations of motion is written as 



R, 



dH 
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dH 
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(13) 



We cool down the system using the following equations 
of motion in which we introduce extra friction terms to 



the above equations [28| : 
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Here, the friction coefficients £r and £p are positive def- 
inite, which determine the relaxation time scale and lead 
to a monotonic decrease of the total energy. 

Instead of the normal kinetic temperature, which loses 
its physical meaning for the system with momentum- 
dependent potentials, we use effective temperature T e g 
proposed by Ref. (29j : 



2 Teff 



AT ^2 1 dt 



(15) 



where J\f is the total number of particles. If we per- 
form Metropolis Monte Carlo simulations with the set- 
ting temperature T sot , the long-time average of the ef- 
fective temperature coincides with T sct quite well [2~ij |. 
This shows T c g is consistent with the temperature in the 
Boltzmann statistics. 

To obtain the equilibrium state at finite temperatures, 
we use the Nose-Hoover thermost at 1 32M-34I] modified for 
momentum-dependent potentials [24l. l35j| . The Hamilto- 
nian of the system with this thermostat is [36[ 
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Here U is the momentum-dependent potential, s is the 
additional dynamical variable for time scaling, p s is the 
momentum conjugate to s, Q is the thermal inertial pa- 
rameter corresponding to a coupling constant between 
the system and the thermostat, g is a parameter to be 
determined as 3iV by the condition for generating the 
canonical ensemble in the classical molecular dynamic 
simulations, and (3 is the reciprocal of T sct of the ther- 
mostat. Then equations of motion are 
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where £ means the thermodynamic friction coefficient. 
During the time evolution described by the above equa- 
tions, Hnosc is conserved and T c ff fluctuates around T sct . 
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III. PURE NEUTRON MATTER AND STABLE 
NUCLEI 

To understand the properties of our QMD models, we 
first calculate the energy and the proton chemical po- 
tential of pure neutron matter at zero temperature. In 
addition, we investigate the surface diffuseness and the 
surface tension of stable nuclei. These are one of the 
key uncertainties that affects the phase diagram at sub- 
nuclear densities To obtain pure neutron 
matter at zero temperature, we use the frictional relax- 
ation method [Eq. (fT4|) ] with the cooling time scale of 
0(1O 3 ) fm/c. We calculate the proton chemical potential 
at zero temperature from the change of the energy by 
inserting a proton into the pure neutron matter. Here 
we relax the position and momentum of the proton with 
fixing those of neutrons (for more details about the pro- 
cedures, see Ref. [23|). To obtain the ground state of 
finite nuclei, we use the conjugate gradient method |37j . 



Energy E n per baryon of pure neutron matter is shown 
in the left panel of Fig. Q] At subnuclear densities, E n of 
both the QMD models exhibits reasonable values com- 
pared with those of other nuclear models. At lower den- 
sities of p n 0.1 fm~ 3 , model 1 gives relatively small 
energy but close to SkM, which gives the lowest energy 
among the other models. At densities below 0.12 fm~ 3 , 
the energy of pure neutron matter for model 2 is larger 
than that for model 1. This tends to prevent neutrons 
from dripping out of nuclei. As we will see later, the 
number density of dripped neutrons for model 2 is in- 
deed smaller than that for model 1. 

According to Ref. (26|, parameter L of the density- 
dependent symmetry energy coefficient also plays an im- 
portant role in determining the density region of the 
pasta phases. This parameter is directly related to the 
derivative of the energy of pure neutron matter with re- 
spect to p n at the normal nuclear density [see Eq. (4) in 
Ref. [H|]: 



L = 3p 



d 

dp n 



(19) 



where e n is the energy density of pure neutron matter. 
The left panel of Fig. [T] shows a larger slope of the energy 
for model 1 than model 2, which leads to a larger value 
of L for model 1. From Eq. ((15)1. we obtain L = 93 MeV 
for model 1 and L = 80 MeV for model 2. This difference 
affects the density at which matter becomes uniform as 
we will discuss in the next section. 

In the right panel of Fig. [TJ we show the proton chem- 
ical potential p p ^ in pure neutron matter at zero tem- 
perature calculated by the QMD models together with 
those by other nuclear models. This result shows QMD 
model 1 gives slightly lower values of p p ^ at high densi- 
ties compared with other nuclear models, whereas model 



2 gives lower values at low densities. As discussed in 

Refs. [13, EH , uncertainty of p p ^ little affects the phase 
diagram of supernova matter because the number density 
of dripped neutrons is very small [39| . 



In Table |TT] we show several quantities related to the 
surface diffuseness and the surface energy of typical heavy 
nuclei, 56 Fe, 90 Zr, 208 Pb, and 238 U, calculated for each 
QMD model. We calculate the surface diffuseness param- 
eter, which, following the spirit of Ref. [4l|, we define as 



Pi 



\dpi/dr\ 



p,n). 



(20) 



Here we have replaced po in the definition of Ref. [4l| . 
which is employed for the semi-infinite system, by the 
central density pi^ n of the finite nucleus. 

We estimate the surface energy a within the framework 
of QMD by subtracting the contributions of the bulk and 
the Coulomb energies from the total binding energy E: 



E - ffcoul - AW(p- m ,Xp) 
AttR 2 



(21) 



Here A and x p are the mass number and the proton frac- 
tion of the nucleus, Seoul is the Coulomb energy of the 
nucleus, and W(p- m ,x p ) is the bulk energy evaluated for 
the central density p- ln of the nucleus. For W(p- m ,x p ), we 
assume the following form, 



W{p- m ,x p ) 



2 \ k 



1 



+4SV x p 



(22) 



where Wv is the binding energy of symmetric nuclear 
matter at po, K$ is the incompressibility, and k- ln = 
(37r 2 p in /2) 1 / 3 and k = (3ir 2 p /2) 1 / 3 are the wave num- 
bers of nucleon at p — p ln and po, respectively. We set 
W v = -16 MeV, K = 280 MeV, S v = 34.6 MeV, and 
po = 0.165 fm -3 for model 1 [H and W v = -16 MeV, 
K = 280 MeV, SV = 33 MeV, and p = 0.168 fm -3 for 
model 2 [H, H3|. The nuclear radius R in Eq. (gTJ is 
defined as 



U EE [ - — 

, 47T p in 



1/3 



(23) 



Table |H] shows that the surface energy a estimated for 
model 2 is systematically higher than that for model 1. 
We also see that the surface diffuseness parameters bi 
of both neutrons and protons are smaller for model 2 
than model 1, which means model 2 yields steeper density 
profile of the nuclear surface. Both of these two facts 
consistently indicate that the nuclear surface energy for 
model 2 is greater than that for model 1. 
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FIG. 1: (Color online) Energy per baryon (left panel) and proton chemical potential fj, p (right panel) in pure neutron matter 
calculated by the QMD models and several other nuclear forces. The line denoted by SLy4 is from Ref. [21(. The other lines 
(FPS, 1', FPS21, and SkM) are from Ref. [H. 



TABLE II: Several quantities of typical heavy nuclei calculated by each QMD model. Eb/A is the empirical binding energy 
from Ref. [40(]. E/A is the binding energy calculated by QMD. pi n is the central nucleon density of the nucleus. b p and b n are 
the surface diffuseness parameter for protons and neutrons, respectively, defined as Eq. (|2U|) . W is the binding energy defined 
as Eq. (|22[) evaluated for the central density and a proton fraction of the nucleus, a is the surface tension defined as Eq. (|21l) . 





56p e 


90 Zr 


i 8 Pb 


238 it 
92 U 


Eg/A (MeV) 


-8.79 


-8.71 


-7.87 


-7.57 


Model 


1 


2 


1 


2 


1 


2 


1 


2 


E/A (MeV) 


-9.09 


-8.77 


-9.25 


-8.68 


-8.66 


-7.65 


-8.31 


-7.38 


Pin (far 3 ) 


0.226 


0.215 


0.213 


0.211 


0.193 


0.163 


0.173 


0.160 


bp (fm) 


4.0 


3.2 


4.2 


3.2 


3.4 


3.1 


3.3 


3.0 


b n (fm) 


4.1 


3.3 


4.0 


3.5 


3.8 


3.4 


3.8 


3.4 


W(p ln ,x p ) (MeV) 


-14.11 


-14.80 


-14.47 


-14.72 


-14.05 


-14.51 


-14.18 


-14.26 


a (MeV fm~ 2 ) 


0.66 


1.13 


1.03 


1.09 


0.62 


1.07 


0.73 


1.03 



IV. NUCLEAR MATTER AT SUBNUCLEAR 
DENSITIES 

A. Simulations and snapshots 

We studied the (n, p, e) system of the proton fraction 
x p = 0.3 at zero and nonzero temperatures. The value 
of x p = 0.3 is typical for matter in collapsing super- 
nova cores. For this purpose, we use 2048 nucleons (614 
protons and 1434 neutrons) in a cubic simulation box 
with periodic boundary condition [HI . We assume that 
the system is not magnetically polarized, i.e., it contains 
equal numbers of protons (and neutrons) with spin up 
and spin down. The relativistic degenerate electrons can 
be well approximated as a uniform background in our sit- 
uations [43-45]. The Coulomb interaction is calculated 
by the Ewald sum (expressions used in our simulations 
are given in Ref. [Hj]). To obtain equilibrium states both 
at zero and nonzero temperatures, we perform simula- 
tions in the following procedure. We first prepare nu- 
clear matter at T — 10 MeV by equilibrating the system 



for about 3000 fm/c using the Nose-Hoover thermostat. 
To reproduce the ground state, we cool down the sys- 
tem with the frictional relaxation (fT4|) for the time scale 
of 0(1O 3-4 ) fm/c. This time scale is sufficiently larger 
than the relaxation time scale of the system O(10 2 ) fm/c, 
which is determined by the typical length scale of the 
structure O(10) fm and the sound velocity O(10 _1 )c. 
To obtain nuclear matter at some fixed nonzero temper- 
ature of Tsot, we start from a snapshot with the effective 
temperature T c ff — T sct obtained in the above cooling 
process. We then equilibrate it with the Nose-Hoover 
thermostat for at least 5000 fm/c. 



In Fig. [21 we show nucleon distributions of the pasta 
phases at zero temperature obtained by the simulations 
for model 2. Compared with those for model 1 (see, e.g., 
Fig. 2 in Ref. (22j), less dripped neutrons are observed. 
Especially, dripped neutrons almost disappear when nu- 
clei become planar. This would be due to the higher en- 
ergy E n of pure neutron matter at low densities as shown 
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FIG. 2: (Color online) Nucleon distribution of the pasta phases at zero temperature for QMD model 2. Simulations are 
performed with 2048 nucleons at a proton fraction x p = 0.3. Each red (blue) particle corresponds to a proton (neutron). Each 
picture shows the pasta phase with (a) spherical nuclei (O.lpo = 0.0168 fm -3 ), (b) cylindrical nuclei (0.2po = 0.0336 fm -3 ), 
(c) slablike nuclei (0.393p = 0.0660fm -3 ), (d) cylindrical holes (0.49p = 0.0823 fm" 3 ), and (e) spherical holes (0.575p = 
0.0966 fm" 3 ). Box sizes are (a) 49.58 fm, (b) 39.35 fm, (c) 31.42 fm, (d) 29.19 fm, and (e) 27.67 fm. 



in Fig. [T] As stated later in Sec. IIV CI this relative lack 
of dripped neutrons decreases the density at which the 
fission instability of spherical nuclei occurs. 



In Fig. [3] we show how slablike nuclei melt to be uni- 
form nuclear matter with increasing temperature. We 
can see a basic picture of phase transitions at nonzero 
temperatures for fixed densities from this figure. Sup- 
pose we increase the temperature from zero. At T ~ 1 
MeV, neutrons start to evaporate from nuclei. At T = 2- 
3 MeV, increase of the volume fraction of nuclear matter 
region due to its thermal expansion triggers a transition 
of the nuclear structure: slablike nuclei start to connect 
with each other to form cylindrical bubbles. At T ~ 4 
MeV, protons cannot be completely confined in nuclei 
and start to evaporate. In this situation, nuclear sur- 
face can no longer be identified, but in many cases there 
remains clustering of protons and neutrons, i.e., nuclear 
matter does not completely become uniform. Finally, the 
clustering inside nuclear matter completely disappears 
and matter becomes uniform. We observe this transition 
at T = 6-7 MeV in the case of Fig. [3] 



B. Identification of phases 

As can be seen from the snapshots of the nucleon dis- 
tribution in the previous section, especially at nonzero 
temperatures, nuclei have complicated shapes and more- 
over nuclear surface becomes diffuse. We thus need to 
quantitatively identify the nuclear shape from nucleon 
distributions obtained by the simulations. For this pur- 
pose, we use the two-point correlation function and the 
Minkowski functionals, especially the area-averaged in- 
tegral mean curvature (H) and the Euler characteristic 
density x/V ■ The Euler characteristic % is a purely topo- 
logical quantity and is expressed as 

X = (number of isolated regions) — (number of tunnels) 
+ (number of cavities). (24) 



For detailed procedures of calculating the Minkowski 
functionals, see Ref. [23[ (see also Refs. [46T - |48j for the 
algorithm of the calculation). 

With these quantities, we can completely classify the 
following typical pasta phases: 

(H) > and x/V > for spherical nuclei (SP) 
(H) > and x/V = for cylindrical nuclei (C) 
< (H) = and x/V = for slablike nuclei (S) (25) 
(H) < and x/V = for cylindrical holes (CH) 
(H) < and x/V > for spherical holes (SH). 

As shown above, x/V is always positive or zero for 
these phases. However, in our previous studies [23l. [24T|. 
and also in the present study, we obtain "spongelike" 
phases with multiply connected structures characterized 
by x/V < 0. We call the phases with x/V < in- 
termediate phases, which appear in the density region 
between those of the phases with rodlike nuclei and slab- 
like nuclei [here we denote as (C,S)], and between slablike 
nuclei and rod-like bubbles [denote as (S,CH)]. The for- 
mer phase (C,S) gives (H) > and x/V < and the 
latter one, (S,CH), gives (H) < and x/V < 0. Consid- 
ering similarity of exotic structures observed in nuclear 
matter and diblock copolymers, several authors pointed 
out a possibility of more complex structures than ordi- 
nary pasta structures, i.e., rods or slabs [IfJ 0, H§| . In 
diblock copolymer melts, one experimentally finds com- 
plicated structures, e.g., so-called gyroid (G) structure 
and the ordered bicontinuous double diamond (OBDD) 
structure. Although the intermediate phases obtained in 
our studies are different from G and OBDD phases, it is 
possible that some complicated structure other than one- 
dimensional lattice of slablike nuclei, hexagonal lattice of 
rodlike ones, and BCC lattice of spherical ones appears 

HE El- 

The quantities (H) and x/V can be calculated only if 
nuclear surface can be identified. Suppose we increase 
the temperature, nuclei start to melt and nuclear sur- 
face cannot be necessarily identified. In this situation, 
the density inhomogeneity of long wavelength starts to 
be smoothed out but still remains. Further increasing 
temperature and exceeding some critical point, matter 
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T=lMeV T=2MeV T=3MeV T=4MeV T=5MeV I=6MeV T=7MeV 




FIG. 3: (Color online) Nucleon distribution at T — 1-7 MeV and a fixed density for model 2. Simulations are performed with 
2048 nucleons at x p = 0.3 and p = 0.393p = 0.0660 fm" 3 (the box size is 31.42 fm), where the phase with slablike nuclei is 
obtained at zero temperature. Each red (blue) particle corresponds to a proton (neutron). 



becomes uniform. To determine the boundary where in- 
homogeneity disappears, we use the two-point correlation 
function defined as 

&( r ) = hj^hj dx ^( x )^( x + r ) ( 26 ) 
= (^(x^x + r))^. (27) 

Here i specifies the species of particles (proton or neu- 
tron, or both collectively), (■ • • } x ,n r denotes an average 
over the position x and the direction of r, and <5j(x) is 
the fluctuation of the density field p^ (x) given by 
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P W(x) 
di{X) = == 



P 



(i) 



(28) 



with the average density of protons or neutrons or both 
collectively: 



v 



(29) 



If the system of nuclear matter is separated into liq- 
uid and gas phases, the two-point correlation function 
£nn(.t) of nucleons oscillates around zero even at long 
distances. Otherwise, £,NN{r) is almost zero (for r > 7 
fm in our cases) without osillating. In this case we judge 
the system is uniform (see Ref. [24| for more details). 



C. Phase diagrams 

In Fig. [4] we show the phase diagram of the pasta 
phases at zero temperature calculated by QMD. The up- 
per and the lower panels of Fig. 3] are the phase diagram 
for model 1 and 2, respectively. The sequence of nuclear 
shapes with increasing density is the same as that pre- 
dicted by all the previous works including those by QMD. 
The density region of the phases with nonspherical nu- 
clei for model 2 is larger than that for model 1: spherical 
nuclei begin to elongate at a lower density and spherical 
bubbles remain until a higher density for model 2. 

The decrease of the density at which nuclei start to be 
deformed would be due to the smaller number density of 
dripped neutrons for model 2 compared with model 1, 
which is originated from the larger energy of pure neu- 
tron matter at low densities shown in Fig. Q] Nuclei are 



SP:sphere CH :cylindrical hole 
C:cylinder SH:spherical hole 
S:slab (, ) intermediate phase 

FIG. 4: Phase diagram of the pasta phases at zero temper- 
ature for model 1 (upper panel) and model 2 (lower panel). 
Proton fraction is x p — 0.3. Horizontal axis is normalized in 
units of the normal nuclear density po for each model. Ab- 
breviations SP, C, S, CH, and SH mean phases with spherical 
nuclei, cylindrical nuclei, slablike nuclei, cylindrical holes, and 
spherical holes, respectively. The parentheses (A,B) show an 
intermediate phase between A and B phases. 



more neutron rich due to the smaller number density of 
dripped neutrons. The nucleon number density inside 
the nuclei is relatively small because of the smaller satu- 
ration density of nuclear matter at the lower proton frac- 
tion. Thus the volume fraction of the nuclei increases, 
which decreases the fission instability of spherical nuclei 
[44| , Furthermore, the decrease of the number density 
of dripped neutrons increases the mass number of nuclei, 
which also increases their volume fraction. As a result, 
although the saturation density of asymmetric nuclear 
matter of x p = 0.3 for model 2 is higher (see below), the 
volume fraction of nuclei for model 2 is about 10 % higher 
than that for model 1 at p = 0.1 pQ. Thus the density 
at which the fission instability occurs for model 2 should 
be ~ 0.0 lpo lower than that for model 1. This value 
is consistent with our present results, which indicate the 
difference between the two models is smaller than 0.05po- 

As to the boundary between the regions of the phase 
with spherical holes and of the uniform phase, however, 
the increase of the density p m at which matter becomes 
uniform would be due to higher saturation density at a 
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proton fraction x p = 0.3 for model 2. The saturation den- 
sity and the saturation energy per baryon at x p = 0.3 are 
0.136 fm -3 and -12.01 MeV for model 1, and 0.147 fm -3 
and —9.86 MeV for model 2 The saturation density 
directly affects p m and, consequently, the p m of model 2 
is higher than that of model 1. Here one should keep in 
mind an effect of the surface energy. In Sec. Mil we show 
that model 2 gives a larger surface tension compared with 
model 1 . The larger surface tension of model 2 favors the 
uniform phase without bubbles and acts to decrease p m . 
However, this effect should be small compared to the con- 
tribution of the saturation density, which can be under- 
stood by taking account of the incompressible property 
of nuclear matter; in the incompressible limit, p m = p s 
and the surface tension does not affect p m at all. 

The result that the model 2 yields a wider density re- 
gion of the pasta phases compared with model 1 can be 
also understood in terms of the density-dependence pa- 
rameter L of the symmetry energy. Within a macroscopic 
model employed in Ref. [26[ [see Eq. (1) and Eq. (4) of 
that reference] , the saturation density at a fixed value of 
x p is given by [26j 

p s {x p )=p [l-3L(l-2x p ) 2 /K ], (30) 

where Kq is the incompressibility. This equation means 
smaller L yields higher p s (x p ) of asymmetric nuclear mat- 
ter. The above higher p s (x p = 0.3) of model 2 than that 
of model 1 shows a smaller L of model 2, which is con- 
sistent with the results of L obtained from the energy of 
pure neutron matter in Sec. IHIl According to a result 
of a comprehensive analysis of Ref. [26j, there is a sys- 
tematic trend that nuclear models with smaller L yields 
a wider density region of the pasta phases. 



Phase diagram of model 2 for x p = 0.3 at nonzero tem- 
peratures is shown in Fig. Each region of this phase 
diagram is defined as follows: (a) SP, (b) C, (c) (C,S), 
(d) S, (e) (S,CH), (f) CH, (g) SH, (B) phase separating 
region, and (C) uniform matter. Abbreviations SP, C, S, 
CH, SH, and (A,B) are the same as in Fig. 2] Compared 
with the phase diagram of model 1 (see Fig. 19 in Ref. 
[24|). both the pasta phases and the liquid-gas phase sep- 
arating region survive until higher temperatures: in the 
case of model 1, the nuclear surface cannot be identified 
above 2-3 MeV, and the critical temperature T c of the 
phase separation is at T > 6 MeV, whereas for model 2, 
the surface melting temperature is at 3-4 MeV and T c 
is at T > 9 MeV. The increase of the surface melting 
temperature can be explained by higher energy of pure 
neutron matter, which prevents neutrons from dripping 
out of nuclei. In addition, as described in Sec. IIII1 the 
surface diffuscness for model 2 is smaller than that for 
model 1 . These properties keep the high density contrast 
between the inside and the outside of nuclei, and hence 
nuclear surface remains to be identified at higher temper- 
atures. However, the increase of the critical temperature 



of the phase separation for model 2 may be explained by 
the smaller value of L, which increases the density where 
the proton clustering instability takes place [2a |. Phase 
separation at high temperatures would be also induced by 
the proton clustering. Thus the same expectation may be 
possible for this situation, i.e., higher symmetry energy 
due to the smaller value of L destabilizes uniform mat- 
ter against the phase separation. As a result, the phase 
separation occurs at relatively higher temperatures for 
model 2. 



V. DISCUSSION AND CONCLUSION 

We investigated the phase diagram of the pasta phases 
both at zero and nonzero temperatures by QMD with two 
different models. Properties of these two models are com- 
pared by calculating the energy and the proton chemical 
potential of pure neutron matter, the surface diffuseness 
and the surface energy of several typical heavy nuclei. 
Differences in the phase diagram, especially the expan- 
sion of the density and temperature region of the nonuni- 
form phases can be explained by these properties of pure 
neutron matter. The sequence of the nuclear shape with 
increasing density and the qualitative feature of thermal 
fluctuation on the nucleoli distribution with increasing 
temperature are the same as observed in our previous 
study [24(. The general picture of the change of the nu- 
cleoli distribution at a fixed density with increasing tem- 
perature is as follows. 

At low temperatures, T — 1-1.5 MeV for model 1 and 
T = 1-2 MeV for model 2; the number density of evap- 
orated neutrons increases with temperature. However, 
the structure of nuclei does not largely change from that 
at T = 0. The nuclear surface becomes diffuse and the 
volume fraction of nuclei increases by thermal expansion. 

At intermediate temperatures, T = 1.5-2.5 MeV for 
model 1 and T — 2-3 MeV for model 2; the nuclear shape 
is significantly deformed and in some cases phase transi- 
tion between different nuclear structures is triggered by 
the increase of the volume fraction of nuclei. Thus the 
density of the phase boundary between the different nu- 
clear shapes decreases with increasing temperature. 

At high temperatures, T ~ 2.5-3 MeV for model 1 
and T ~ 3-4 MeV for model 2; evaporated nucleons are 
dominant and the nuclear surface can no longer be iden- 
tified. However, the long-range correlations between nu- 
cleons due to the liquid-gas phase separation remain until 
a higher temperature T ~ 6 MeV for model 1 and T ~ 9 
MeV for model 2. Above these temperatures, inhomo- 
geneity disappears at any density. 

The density-dependence parameter L of the symme- 
try energy is the key to understand the uncertainties 
of the density region of the pasta phases in cold neu- 
tron stars [2y|. This parameter is also helpful to un- 
derstand the present results and to predict the general 
tendency of phase diagram of the pasta phases in super- 
nova cores. From the energy of cold neutron matter, we 
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FIG. 5: (Color online) Phase diagram of nuclear matter of x p = 0.3 at subnucear densities by QMD model 2 plotted in the 
p-T plane. The horizontal axis is normalized in unit of the nuclear saturation density. The dashed lines correspond to phase 
separation lines. The dotted lines show the boundary above which nuclear surface cannot be identified. The dash-dotted lines 
show boundaries between different phases. Abbreviations are the same as described in the caption to Fig. [4] The meanings of 
regions (a)-(g), (B), and (C) are explained in the text. Simulations have been carried out at the points denoted by circles. 



obtain L = 93 MeV for QMD model 1 and L = 80 MeV 
for model 2. The larger critical temperature T c for the 
phase separation of model 2 would be due to the smaller 
value of L. In addition, the smaller L also yields higher 
saturation density of asymmetric nuclear matter [see Eq. 
(|30[) ], which in turn increases the density at which the 
system becomes uniform nuclear matter. If one uses a 
nuclear force with a smaller value of L, the density and 
temperature region of nonuniform nuclear matter would 
broaden 

Let us now discuss the abundance of nonuniform 
phases of nuclear matter in supernova cores. Here we 
take EOS by Shen et al. [5l[ as an example, which is 
one of the widely used EOSs in supernova simulations. 
For the nuclear model employed in this EOS, we estimate 
L ~ 120 MeV, which is rather high in the range of un- 
certainty of nuclear forces (see, e.g., Fig. 1 of Ref. [26}). 
In our previous work Q, using EOS by Shen et al, we 



have estimated the mass fraction of the pasta phases just 
before bounce and have obtained ~ 10-20 %. Because we 
can expect values of L for other typical EOSs are smaller 
than L — 120 MeV for EOS by Shen et al, our previous 
estimate of the mass fraction using Shen's EOS would be 
very conservative. It is reasonable to conclude that the 
mass fraction of the pasta phases would be larger than 
10-20 % in the later stage of the collapse. 

As we have shown in Ref. [8j, neutrino opacity via 
weak neutral current in the pasta phases can be signifi- 
cantly different from that without taking account of the 
pasta phases. Furthermore, even if nuclear surface melts, 
the neutrino opacity of the vector current contribution in 
the liquid-gas phase separating system is still larger than 
that of the completely uniform gas phase. Our present 
result also indicates expansion of the phase separating 
region for smaller values of L. Although the influence 
of the phase separating region on the mechanism of the 
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supernova explosion has yet to be revealed completely, it 
could exist not only during collapsing phase but also after 
bounce. Depending on the location of the neutrino pho- 
tosphere and the temperature profile in the late stage, 
these region can affect the success of supernova explo- 
sions 153 1. 
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